% fiven initial condition, total flight time and the normal vector of the
% Sun, calculate the visibie and invisible region of the Sun.
% IC = initial conditions
% T  = period of the periodic orbit
% FN = normal vector lists.


function plotVisibility(IC,T,u,FN,CameraAngleCapability)
options=odeset('RelTol',1e-13,'AbsTol',1e-13);

to = [0,T];
Xo = IC;
[t,S] = ode113(@(t,S)CR3BP_n(t,S,u),to,Xo,options);
S = S';
% plot3(S(1,:),S(2,:),S(3,:))
[row,col,~] = size(FN);
visibility=zeros(row,col);

parfor ii = 1:row
    for jj = 1:col
        normvec = [FN(ii,jj,1);FN(ii,jj,2);FN(ii,jj,3)];
        for k = 1:length(t)
            posvec = S(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if angle_temp < CameraAngleCapability
                visibility(ii,jj) = 1;
                break
            end
        end
    end
end
ImageBinary = mat2gray(visibility);
% figure;
% imshow(ImageBinary)
Window = figure('OuterPosition',[100,100,1500,800]); hold on;
mesh(ImageBinary,'EdgeColor','none','FaceColor','flat');
xlabel('Longitude')
ylabel('Latitude')
title('Sun Visibility')
xticks(linspace(0,row,5));
xticklabels({'180W','90W','0','90E','180E'});
yticks(linspace(0,col,5))
yticklabels({'90S','45S','0','45N','90N'});
xlim([1,row]);
ylim([1,col])

Window = figure('OuterPosition',[100,100,1500,800]); hold on;
contour(ImageBinary,[1,1]);
xlabel('Longitude')
ylabel('Latitude')
title('Sun Visibility')
xticks(linspace(0,row,5));
xticklabels({'180W','90W','0','90E','180E'});
yticks(linspace(0,col,5))
yticklabels({'90S','45S','0','45N','90N'});
xlim([1,row]);
ylim([1,col])
end
